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Abstract 

Developments in dynamical systems theory provides new support 
for the discretisation of pdes and other microscale systems. Here 
we explore the methodology applied to the gap-tooth scheme in the 
equation-free approach of Kevrekidis in two spatial dimensions. The 
algebraic detail is enormous so we detail computer algebra procedures 
to handle the enormity. However, modelling the dynamics on 2D 
spatial patches appears to require a mixed numerical and algebraic 
approach that is detailed in this report. Being based upon the com- 
putation of residuals, the procedures here may be simply adapted to 
a wide class of reaction-diffusion equations. 
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1 Introduction 

We extend the dynamical systems, holistic, approach to the macroscale dis- 
crete modelling [8, 10, 11, e.g.] to two dimensional, homogeneous, non- 
linear reaction-diffusion equations on coupled patches of space. Following 
the 'equation free' approach of Kevrekidis and colleagues [5, e.g.], we ad- 
dress the extraction, using dynamical systems theory, of computationally 
efficient macroscale models from given microscopic models, whether PDE or 
lattice dynamics or other microscale systems. Here we bridge space scales by 
generalising to multiple dimensions (specifically 2D) the methodology and 
supporting theory for the 'equation free', gap-tooth method for microsimu- 
lators [3, 13, 14]. As a particular example, this report uses the real valued, 
two dimensional, Ginzburg-Landau equation 

_ = V 2 u+ ct(u-u 3 ) . (1) 
ot 

We choose this 2D real Ginzburg-Landau equation as a prototype PDE be- 
cause it is well studied and its dynamics well understood [4, 7, e.g.]. This 
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report details the construction by computer algebra of the macroscale dis- 
crete model of its dynamics in two spatial dimensions. The general theo- 
retical support, the performance and the physical interpretation are detailed 
elsewhere. 

Place the discrete modelling of two dimensional, reaction-diffusion equa- 
tions within the purview of centre manifold theory by dividing the domain 
into square patches and introducing special interpatch coupling conditions. 
Define a grid of points (Xi,Yj) with, for simplicity, constant spacing H. The 
i, jth patch, Ey, is then centred upon (Xi, Yj) and of width Ax = Ay = 2rH: 
when the ratio r < 1/2, Ey forms a patch separated by empty space (gaps) 
from neighbouring patches; when the ratio r = 1 the patches overlap and 
the analysis reduces to that of holistic discretisation [8, 10, 11, e.g.]. Define 
that Uy (x, y, t) denotes the field in the i, jth element and so evolves accord- 
ing to the Ginzburg-Landau pde (1). Using the parameter y to control the 
strength of the coupling, use coupling conditions around the i, jth patch of 

uy (Xi ± tH, i, , t) = St r ^Kj (Xi, Yj , t) , |y - Y, I < rH , 

Uy(x,Yj ±rH,t) = £^(y)£± r (y)uy (X^t) , |x-X t | < rH, 

In the overlapping element case, r = 1 , we use 'interelement coupling condi- 
tions' around the i, jth element of 

Uy(X i±1 ,-y,t) =yUi±y(X i±1 ,-y,t) + (1 -Y)uy(Xi,-y,t), |y - Y } | < H, 
u-y(x, Y j±1 ,t) = yuy ±1 (x, Y j±1 , t) + (1 -y)uy(x,Yj,t), |x-X t | < H. 

(3) 

These are a natural extension to 2D of those established for ID dynamics [12]. 
Then centre manifold theory [1, 2, 6, e.g.] assures us of the existence, rele- 
vance and approximation of a slow manifold, macroscale model parametrised 
by a measure of the amplitude in each element and the coupling strength y. 

However, it appears that models of the pde (1) cannot be constructed 
algebraically. Thus Section 4 details computer algebra to numerically solve 
for the microscopic subgrid scale field. That is, we actually create a coarse 
grid model of a the fine grid/lattice dynamics of an approximation to the pde. 
Its application here to the discretisation of the Ginzburg-Landau pde (1) 
serves as a proof of principle for applying holistic discretisation to general 
PDEs of two or more spatial dimensions. 

For example, consider the pde (1) when discretised inside patches by a 
5x5 microscale lattice (n = 2). Using the patch ibcs (2) the slow manifold 
discretisation is 

Uy = ^6 2 Uy - t^O - t)*Xj + ^ (l - £) (l - S) 
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+ a(Uy-U y ) 

+ KY^Ys { [ 3( ^ 5xUi 'i )2 + K 5 xUy) 2 ] (2 + 6 2 ) + (5 2 Uy) 2 

+ [3(^5,Uy) 2 + l(S 2 Uy) 2 ] (i + s^ + ^Uyfliiy 
+ (9(a 2 + y 4 ), (4) 

where the bold centred difference operator applies in both spatial dimensions, 

6 2 Uy = Ui+u+Ui-u+Uy+T+Uy-T^Uy, (5) 
6 4 Uy = U i+2,j + Ui-2,j + Uy +2 + Uy_ 2 

- 4(U i+1>j + Ui_y + Uy +1 + Uy_! ) + 12Uy , (6) 

2 Construct the algebraic slow manifold 

This section provides computer algebra code to generate analytic holistic 
discretisations of reaction-diffusion pdes in 2D. Different pdes are analysed 
by changing the nonlinear term. Higher order models are constructed by 
changing the order of neglected terms. Unfortunately, analytic construction 
can only be carried out to low order accuracy: Section 4 describes code for a 
numerical construction that also applies to patches in the gap-tooth scheme. 
All code is written in the computer algebra package Reduce. 1 

2.1 Initialisation 

Set some parameters to improve printing of the results. 

t» casm2d « 

% see cadsmpd2d.pdf for documentation 
on div; off allfac; on revpri; 
factor gam,h,eps,nu,alf ; 

o o o 

Subgrid variables The subgrid, intra-element, structures are functions of 
intra-element microscale variables xi = (x — X{)/h and yi = [y — yj)/H. 

» casm2d «+ 

depend xi , x ; 

let df (xi,x)=>l/h; 

depend yi,y; 

^ttp : / / www . reduce- algebra . com 
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let df (yi,y)=>l/h; 

o o o 

Parametrise the slow manifold Parametrise the slow manifold, macroscale, 
evolution by the evolving grid values u(i , j ) = such that dlly/dt = gi j . 

» casm2d «+ 

operator u; 
depend u,t; 

let df (u(~k, ~m) ,t)=>sub(-[i=k, j=m} ,gij ) ; 

o o o 

The linear, slow subspace, approximation is that of piecewise constant 
fields and no evolution. 

» casm2d «+ 

uij :=u(i,j) ; 
gij:=0; 

o o o 

Set asymptotic truncation Here scale nonlinearity parameter a with 
parameter y so we truncate to residuals and errors of order (9 (a 3 +y 3 ). 

» casm2d «+ 

let eps~3=>0; 
gamma : =gam*eps ; 
alpha : =alf *eps ; 

o o o 

General multinomial For the method of undetermined coefficients we set 
up a general multinomial solution with unknown coefficients cc(m,n), up to 
order o in the intra-element variables. Collect the unknown coefficients in 
the set cs. 

» casm2d «+ 

o : =6 ; 

operator cc; 

vv:=for m:=0:o sum for n:=0:o-m sum cc(m,n)*xi~m*yi~n$ 
cs:=for m:=0:o join for n:=0:o-m collect cc(m,n)$ 

o o o 

Need to find the coefficients in multivariate expressions so define a recur- 
sive procedure coef f s that given a list of expressions and a list of variables 
it creates a new extended list of all the coefficients of the variables. 
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» casm2d «+ 

procedure coef f s (exps , vars) ; 

if vars={} then exps else coeffs( 

foreach z in exps join coef f(z, first vars), 
rest vars) ; 

o o o 

2.2 Iteratively construct the slow manifold 

Now start the iteration, repeating corrections until all residuals are zero (to 
a maximum of ten iterations). 

» casm2d «+ 
for it: =1:10 do begin 

o o o 

Compute residuals of governing equations Compute the residual of 
the reaction-diffusion pde in the general i, jth element. Could easily change 
the reaction term to any polynomial in the field uy . Could also incorporate 
advection terms into the microscale pde. Also compute the residuals of the 
inter-element coupling ibcss. 

t» eqnResiduals « 

de : =df (uij , t)-df (uij , x, 2) -df (uij ,y, 2) -alpha* (uij-ui j "3) $ 
bcr :=sub(xi=l,uij)-sub(xi=0,uij) 

-gamma* (sub ({xi=0 , i=i+l} , uij ) -sub (xi=0 , uij ) ) $ 
bcl : =sub(xi=-l ,uij )-sub(xi=0,uij ) 

-gamma* (sub ({xi=0 , i=i-l} , uij ) -sub (xi=0 ,ui j ) ) $ 
bet : =sub (yi=l , uij ) -sub (yi=0 , uij ) 

-gamma* (sub ({yi=0 , j = j +1} , uij ) -sub (yi=0 , uij ) ) $ 
beb : =sub (yi=-l , uij ) -sub (yi=0 , uij ) 

-gamma* (sub ({yi=0 , j=j -1} , uij ) -sub (yi=0 , uij ) ) $ 

ooo 

Include the computation of the residuals both here and later. 

» casm2d «+ 

« eqnResiduals t» 

For information as to the progress of the iteration, print out the lengths 
of the residuals: when all are one, then all residuals are probably zero. 
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» casm2d «+ 

write lengths :={ length (de) , length (bcr) , length (bcl) 
, length (bet) , length (beb)}; 

o o o 

Add the as yet unknown corrections To find the desired update in each 
iteration, first substitute the form of the update into the computed residuals. 

» casm2d «+ 

deq : =de+gd-df (vv , x , 2) -df (vv , y , 2) ; 
rbcr : =bcr+sub (xi=l , vv) -sub (xi=0 , vv) ; 
lbcl : =bcl+sub (xi=-l , vv) -sub (xi=0 , vv) ; 
tbet : =bct+sub (yi=l , vv) -sub (yi=0 , vv) ; 
bbcb:=bcb+sub(yi=-l ,vv)-sub(yi=0,vv) ; 

o o o 

Solve for the corrections Then extract equate coefficients of each power 
of the multinomial in the intra-element variables, and append boundary con- 
ditions. 

» casm2d «+ 

eqns:=append(coeffs({deq},-[xi,yi}) ,cc(0,0) . 

append (coeff (rbcr ,yi) , append (coeff (lbcl ,yi) , 

append ( coef f (tbet ,xi) , coeff (bbcb,xi))))) ; 
sol : =solve (eqns , gd . cs) ; 

o o o 

Update Update the field and the evolution, assuming a solution was found 
(not true for higher orders). 

» casm2d «+ 

uij : =uij+sub(sol , vv) ; 
gij :=gij+sub(sol,gd) ; 

o o o 

Terminate End the iteration when all residuals are zero, or too many 
iterations have been performed. 

» casm2d <M+ 

showtime ; 

if {de,bcr,bcl,bct,bcb}=-[0,0,0,0,0} then write it :=1000000+it; 
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end; 



o o o 



2.3 Scrounge an extra order of evolution using solv- 
ability 

First define the linear integral operator inthat (a, xi) = — |£j)ad£, in 
order to quickly apply the solvability condition. 

t» casm2d «+ 

operator inthat; 
linear inthat; 

let { inthat (~a~~p, ~b)=>0 when (a=b)and(not evenp(p)) 
, inthat (~a~~p,~b)=>2/(p+l)/(p+2) when (a=b)and evenp(p) 
, inthat (~a, ~b)=>0 when (a=b) 
, inthat (l,~b)=>l 

>; 

o o o 

Set the requisite next order of truncation in the asymptotic expansion by 
finding the highest order currently retained in coupling y, then setting to 
discard terms of two orders higher. Use the 'instant evaluation' property of 
for-loop variables in Reduce. 

» casm2d «+ 

o : =2+deg( (1+eps) "9 , eps) $ 
for p:=o:o do let eps~p=>0; 

o o o 

Compute exactly the same residuals of the pde and the inter-element 
coupling conditions. 

t» casm2d «+ 

« eqnResiduals t» 

o o o 

Compute the next correction to the evolution gij by integrating the 
residual of the pde over an element, and including the contributions from 
the boundary coupling residuals. 



gd 



t» casm2d «+ 

=inthat(de,xi)$ 
=inthat(gd,yi)$ 

=gd+inthat (bcr+bcl ,yi) /h~2+inthat (bct+bcb ,xi) /h"2 ; 
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gij :=gij-gd$ 
showtime ; 

ooo 



3 Postprocessing analytic 

3.1 Write a sci/matlab function script 

Optionally generate an efficient sci/matlab function in file gl2dnumred. This 
function file is to be processed with the unix script red2mat and then it can 
be used in matlab to integrate solutions using ode solvers. 

» casm2d «+ 

write " 

Generating mat/scilab function file in gl2dnumred 

ii . 

> 

« functionfile » 

write "finished generating mat/scilab file"; 

write "Use sed -f red2mat gl2dnumred > gl2dnpde.m"; 

showtime ; 

ooo 

Improve printing. » functionfile « linelength 66$ off nat$ ooo 

Need boundary conditions that depend upon the width of the stencil of 
the discretisation. Get this width from the truncation in parameter eps, 
assuming coupling y is proportional to eps. If not, then have to change. 
» functionfile «+ o:=deg((l+eps)~9,eps) ; ooo 

Want code for the fully coupled system, y = 1 , and without the artificial 
parameter eps. » functionfile «+ gam:=eps:=l; ooo 

Use the scope package to generate efficient computation of the discreti- 
sation. » functionfile «+ load_package scope; ooo 

Now direct output to the file gl2dnumred. » functionfile «+ out 
"gl2dnumred"$ ooo 

Write the preamble for the function. Given a vector of macroscale vari- 
ables, unpack them into an array of 2D grid values. Get the grid size and 
the gap-tooth ratio as global variables. 

» functionfile «+ 

write "function udot=gl2dnpde(t ,uv) ; 
•/,'/, Tony Roberts, Jan 2011 
global h rr; 
nn=length(uv) ; 
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n=sqrt (nn) ; 

u=zeros (n,n) ; u( : )=uv(l :nn) ; "$ 

o o o 

Allow for two possible boundary conditions: either doubly periodic; or 
odd-even solutions. Choose in the reduce code by changing the switch in this 
if-statement. 

t» f unctionf ile «+ 

if then % use periodic BCs or odd-even BCs 
« periodicBCs »else « oddevenBCs >t>$ 

o o o 

Pad the data array with the periodic BCs: the two cases depend upon 
the stencil width. 

» periodicBCs « 

if o=l then write "°/„ periodic BCs 
u=[u(end, :) ;u;u(l, :)] ; 
u=[u(: ,end) ,u,u(: ,1)] " 
else write "°/ periodic BCs 
u=[u(end- M ,o-l, " :end, : ) ;u;u(l : " ,o, " , : )] ; 
u=[u(: ,end- M ,o-l, " :end) ,u,u( : , 1 : " ,o, ")] " 

o o o 

Alternatively, pad the data array with the odd/even BCs: the two cases 
depend upon the stencil width. The odd/even BCs suit a sine expansion in x 
and a cosine expansion in y. 

» oddevenBCs « 

if o=l then write "°/„ odd-even, sine-cosine BCs 
u=[-u(l, : ) ;u;-u(end, : )] ; 
u= [+u(: ,1) ,u,+u(: ,end)] 11 
else write "°/ odd-even, sine-cosine BCs 
u=[-u(",o,":-l:l, :);u;-u(end:-l:end-",o-l, M , :)] ; 
u= [+u(: ,",o,":-l:l) ,u,+u(: , end : -1 : end-" , o-l , " )] " 

ooo 

Define index vectors that apply to every macroscale grid point in the inte- 
rior of the padded data array. » functionfile «+ write " j = " ,o, " + (1 :n) ; 
i=j"$ ooo 

Use the scope package optimisation to generate the code for the com- 
putation of an array of time derivatives. » functionfile <M+ optimize 
udot: = :gij iname o$ ooo 
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Write code to form the array into a vector. » f unctionf ile «+ write 
"udot=udot(:)"$ ooo 

Close the function file as all code is done. » f unctionf ile «+ shut 
"gl2dnumred" ; ooo 

Undo a printing option, and clear gam and eps so we can see the coupling 
parameter influence again. » f unctionf ile <W+ on nat; clear gam; clear 
eps; ooo 

3.2 Derive the equivalent PDE 

Optionally derive the equivalent PDE for this model. 

» casm2d «+ 
if 1 then begin write " 

Optionally find the equivalent partial differential equation 

of the holistic discretisation. 
11 . 

« equivpde » 
showtime ; 

end; 

ooo 

Depending upon the order of inter-element coupling, as measured by eps, 
then set the order of truncation in the grid size h. This code assumes gamma 
is first order in eps, and needs changing if otherwise. Use a dummy for- loop 
because, in Reduce, loop parameters, here q, have the "instant evaluation" 
property that means it can be used usefully as an index in the let statement. 

» equivpde « 

o:=2+2*deg((l+eps)~9,eps) ; 
for q:=o:o do let h~(q+l)=>0; 

ooo 

Use operator uu(p,q) to denote the mixed derivative 3 p+q U/3x p 3ij C| . 
Hence define a rule to expand u(i+k } in a multivariable Taylor se- 

ries. Reduce does not like 0~0 so handle the zero cases separately in the 
summation. 

» equivpde «+ 

operator uu; 

rules:={u(~k,~D=> uu(0,0) 

+(for q:=l:o sum uu(0 , q) * (1-j ) ~q*h~q/f actorial (q) ) 
+(for p:=l:o sum uu(p,0)*(k-i) ~p*h~p/f actorial (p)) 

11 



+(for p:=l:o sum for q:=l:o-p sum 

uu(p , q) * (k-i) ~p* (1- j ) ~q*h~ (p+q) /factorial (p) /factorial (q) ) 

}$ 

o o o 

Make different printing. » equivpde «+ off revpri; on list; 

o o o 

Transform the discrete model to its equivalent PDEusing the Taylor ex- 
pansion. Remove the artificial expansion parameter eps by settng to one. 
Sum at full coupling gam; remove this to see the y dependence. 

» equivpde «+ 

write 

gpde:=(sub({eps=l,gam=l},gij) where rules) $ 
off list; 

ooo 



3.3 Derive the centred finite difference form 

Optionally derive the finite difference form of the evolution. 

» casm2d «+ 
if 1 then begin write " 

Generate evolution in terms of centred finite difference 
operators md=mu*delta and dd=delta~2 in each of the x and y 
directions . 
> 

« finitediff » 
showtime ; 

end; 

ooo 

Define the four centred difference operators, two in each spatial direction. 
Then define commutation rules so that expressions involving these operators 
change to a canonical form where x operators come before y operators, and 
p.5 operators comes before S 2 operators. Lastly, invoke the operator rule that 
y} = 1 + 8 2 /4. 

» finitediff « 

operator mdx; operator mdy; 
operator ddx; operator ddy; 
linear mdx; linear mdy; 
linear ddx; linear ddy; 
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let { mdy (mdx (' 


'a, 


~u), 


~v) 


=>mdx(mdy (a, v) ,u) 


, ddx (mdx ( ' 


'a, 


~u), 


~v) 


=>mdx(ddx(a, v) ,u) 


, ddy (mdx ( ' 


'a, 


~u), 


~v) 


=>mdx(ddy (a, v) ,u) 


, ddx (mdy ( ' 


'a, 


~u), 


~v) 


=>mdy (ddx (a, v) ,u) 


, ddy (mdy ( ' 


'a, 


~u), 


~v) 


=>mdy (ddy (a, v) ,u) 


, ddy (ddx (' 


'a, 


~u), 


~v) 


=>ddx(ddy (a, v) ,u) 


, mdx (mdx ( ' 


'a, 


~u), 


~u) 


=>ddx(a,u)+ddx(ddx(a,u) ,u)/4 


, mdy (mdy ( ' 


'a, 


~u), 


~u) 


=>ddy(a,u)+ddy(ddy(a,u) ,u)/4 



>; 



o o o 

Define the shift rules that Ui±i = ± (x5Ut + ^5 2 Ui. Apply them to the 
four cases of 11^^ and U i( j ±k . 

» finitediff «+ 

rules : =-[u(i , j )=>u 

,u(i+~k,~j)=>u(i+k-l , j)+mdx(u(i+k-l, j) ,u) 

+l/2*ddx(u(i+k-l, j) ,u) when k>0 
, u (i+~k , ~ j ) =>u (i+k+1 , j ) -mdx (u ( i+k+1 , j ) , u) 

+l/2*ddx(u(i+k+l, j) ,u) when k<0 
,u(~i, j+~k)=>u(i, j+k-l)+mdy(u(i, j+k-1) ,u) 

+l/2*ddy(u(i, j+k-1) ,u) when k>0 
,u(~i, j+~k)=>u(i, j+k+l)-mdy(u(i, j+k+1) ,u) 

+l/2*ddy(u(i, j+k+1) ,u) when k<0 

}; 

o o o 

Transform the evolution of the grid values. Could also transform the 
microscale structure in Uy, but have not done so here. 

» finitediff «+ 

write gop : =(sub(eps=l ,gij ) where rules); 

o o o 

3.4 Fin analytic 

Give the overall grand final 'end' statement, then output all the code. 

» casm2d «+ 

end; 

o o o 
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4 Numerically construct the slow manifold 



This section describes computer algebra code to generates macroscale models 
of the gap-tooth scheme applied to react ion- diffusion pdes in 2D. Different 
pdes are analysed by simply changing the nonlinear term, or more carefully 
changing the Laplacian. Higher order models are constructed by changing 
the order of neglected terms. Instead of trying to solve analytically for the 
subgrid scale structure, here solve for the structure numerically. 

All code is written in the freely available computer algebra package Re- 
duce. 2 

4.1 Initialisation 

Improve the format of printing. 

» ncsm2d « 

% see cadsmpd2d.pdf for documentation 
on div; off allfac; on revpri; 
factor gam,h,eps,nu,alf ; 

o o o 

Load routines to do the LU decomposition, Section 6.1, and subsequent 
back substitutions, Section 6.2. 

» ncsm2d «+ 

in "lu_decomp.red"$ 
in "lu_backsub . red"$ 

o o o 

Need to find the coefficients in multivariate expressions so define a recur- 
sive procedure coef f s that given a list of expressions and a list of variables 
it creates a new extended list of all the coefficients of the variables. 

» ncsm2d «+ 

procedure coef f s (exps , vars) ; 

if vars={} then exps else coeffs( 

foreach z in exps join coef f(z, first vars), 
rest vars) ; 

o o o 

Use patches of size ratio r where: r = 1 is holistic overlapping; r = 1/2 
the elements abut; and for r < 1 /2 there is a gap in the gap-tooth scheme. 



2 http : / / www . reduce- algebra . com 
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» ncsm2d «+ 

%%r:=l; % the relative size of a patch or element 

ooo 

The subgrid, microscale, numerical resolution is improved by increas- 
ing n. The subgrid grid step is dx = H/n and the number of equations and 
unknowns is nn = N = (2n + 1 ) 2 + 1 . Takes some minutes to construct 
solutions for lattice n = 8, but it is reasonably doable; have not tried n > 9 
yet. 

» ncsm2d «+ 

in: =2; % number of dynamic points per half side 
dx:=r*h/n$ % microgrid step size, scaled by factor of r 
ns:=2*n+l$ °/„ number of variables per side of 2D patch 
nn:=ns~2+l$ % total number of equations and variables 

ooo 

Define the matrices used to store the subgrid, and to represent and solve 
the equations. The scope of matrices are global in Reduce. 

» ncsm2d «+ 
matrix eqns(nn, 1) , zeros (nn, 1) ; 

matrix indx(nn, 1) ,vv(nn, 1) ,lu(nn,nn) ; % used by LU 

ooo 

Parametrise the slow manifold The slow manifold of the macroscale 
discretisation is to be parametrised by the evolution of Uy = u(i,j). The 
evolution dlly/dt = gy = gij and unknown updates are stored in gd. 

» ncsm2d «+ 

operator u; depend u,t; 

let df (u(~k, ~m) ,t)=>sub(-[i=k, j=m} ,gij ) ; 

gij :=gd/dx~2; 

ooo 

Initialise with constant field in uij = uy. Also add in the unknown 
update field ud. Use ii and jj as the subgrid microscale lattice indices. 

» ncsm2d «+ 

matrix uij(ns,ns)$ 
operator ud; 

for ii:=l:ns do for jj:=l:ns do 
uij (ii, j j) :=u(i, j)+ud(ii, j j)$ 

ooo 
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Possibly include a forcing function f u (p , q) =df (f f , x , p , y , q) , the phys- 
ical derivatives, evaluated at the (i, j) grid point. But omit for the moment. 

» ncsm2d «+ 

operator fu; 
maxf : =0 ; 

ffij:=(for m:=0:maxf sum h~m 

*(for n:=0:m sum fu(n,m-n)*xi~n*yi~(m-n))) ; 

ooo 



Set asymptotic truncation Truncate the asymptotic expansion in pow- 
ers of the coupling parameter y = gam. Here scale the nonlinearity parameter 
ct with y to most easily control the truncation. 

» ncsm2d «+ 

let { eps~4=>0 }; 
gamma : =gam*eps ; 
gammad : =l-gamma$ 

alpha :=alf*eps~2; % parametrises strength of nonlinearity 
nu:=l; % dissipation parameter 

ooo 

Also ignore the 'updates' when multiplied by the small parameter so that 
the equations remain linear in the updates. 

» ncsm2d «+ 

let { eps*ud(~i , ~j)=>0, eps*gd=>0, 

gam*ud(~i , ~ j )=>0, gam*gd=>0 }; 

ooo 



4.2 Classic Lagrange interpolation 

Coupling conditions are written in terms of classic Lagrange interpolation. 
Thus set up the interpolation fields here as functions of the surrounding 
macroscale grid values. As written so far we can go to errors C(y 5 ). Suspect 
I should be able to code this better. 

First define procedures for centred mean and difference operators in each 
of the two spatial directions. 

» ncsm2d «+ 

procedure mudi(u); (sub(i=i+l ,u)-sub(i=i-l ,u) )/2 ; 
procedure mudj (u) ; (sub(j=j+l,u)-sub(j=j-l,u))/2; 
procedure dedi(u); (sub(i=i+l,u)-2*u+sub(i=i-l,u)) ; 
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procedure dedj (u) ; (sub(j=j+l,u)-2*u+sub(j=j-l,u)) ; 

o o o 

Second, interpolate in the x direction where xi = (x — Xy )/H. 

» ncsm2d <M+ 

uu:=u(i, j)$ 

uu : =uu+gamma* (+mudi (uu) +dedi (uu) *xi/2) *xi 

+gamma"2*dedi (+mudi (uu) +dedi (uu) *xi/4) *xi* (xi "2-1) /6 
+gamma"3*dedi (dedi (+mudi (uu) +dedi (uu) *xi/6) ) 
*xi*(xi~2-l)*(xi"2-4)/120 

+gamma~4*dedi (dedi (dedi (+mudi (uu) +dedi (uu) *xi/8) ) ) 

*xi* (xi"2-l) * (xi"2-4) * (xi"2-9) /5040 

$ 

o o o 

Third, interpolate in the y direction where yi — [y — Yy)/H. 

» ncsm2d «+ 

uu : =uu+gamma* (+mudj (uu) +dedj (uu) *yi/2) *yi 

+gamma~2*dedj (+mudj (uu)+dedj (uu)*yi/4)*yi*(yi~2-l)/6 
+gamma"3*dedj (dedj (+mudj (uu)+dedj (uu)*yi/6)) 
*yi*(yi~2-l)*(yi"2-4)/120 

+gamma"4*dedj (dedj (dedj (+mudj (uu)+dedj (uu)*yi/8))) 

*yi* (yi"2-l) * (yi"2-4) * (yi~2-9) /5040 

$ 

o o o 

4.3 Iteratively construct the slow manifold 

Repetitively find updates to the subgrid microscale field until all residuals 
are zero. 

» ncsm2d «+ 
for iter: =1:19 do begin 

o o o 

Update the u field and evolution Form the residuals of lattice equations 
and all the coupling conditions, using q to count the number of equations. 
First the corner equations. 

» ncsm2d «+ 
write "Compute the u residual"; 
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q:=0; °/ counts the number of equations used 



eqns (q: =q+l , 1) 
eqns (q: =q+l , 1) 
eqns (q: =q+l , 1) 
eqns (q: =q+l , 1) 



=ud(l,l); 
=ud(l ,ns) ; 
=ud(ns , 1) ; 
=ud(ns ,ns) ; 

ooo 



Now adjoin the residuals of the interpatch/element coupling conditions. 
First, if r = 1, then use the interelement coupling conditions. 

» ncsm2d «+ 

if r=l then for ll:=2:ns-l do begin 

eqns(q:=q+l , 1) :=uij (ns , 11)- (1-gamma) *uij (n+1 , 11) 

-gamma*sub(i=i+l ,uij (n+1 , 11) ) ; 
eqns(q:=q+l , 1) :=uij (1 , 11) - (1-gamma) *uij (n+1, 11) 

-gamma*sub(i=i-l ,uij (n+1 , 11) ) ; 
eqns(q:=q+l , 1) :=uij (11 ,ns)- (1-gamma) *uij (11 , n+1) 

-gamma*sub(j=j+l ,uij (11 ,n+l) ) ; 
eqns(q:=q+l , 1) :=uij (11 , 1) - (1-gamma) *uij (11, n+1) 

-gamma*sub(j=j-l ,uij (11 ,n+l) ) ; 

end 

ooo 

For the following patch coupling conditions (when r < 1), recall that 
uu stores an expression for classic Lagrange interpolation from neighbouring 
grid values. 

» ncsm2d «+ 

else for ll:=2:ns-l do begin 

X y :=r*(ll-n-l)/n; % fraction distance along side 



eqns(q: =q+l , 1) 
eqns(q:=q+l , 1) 
eqns(q:=q+l , 1) 
eqns(q:=q+l , 1) 
end; 



=ui j (ns , 11) -sub ({xi=r , yi=xy} , uu) 
=ui j (1,11) -sub ({xi=-r , yi=xy} ,uu) 
=uij (11 ,ns)-sub({yi=r ,xi=xy} ,uu) 
=uij (11, l)-sub({yi=-r,xi=xy}-,uu) 



ooo 



Adjoin the residuals of the interior subgrid equations which are a simple, 
second order, centred difference approximations to the pde. 

» ncsm2d «+ 

for ii:=2:ns-l do for jj:=2:ns-l do 
eqns(q:=q+l,l) :=-df (uij (ii, j j) ,t)*dx~2 

+nu*((uij (ii+1, j j)-2*uij (ii, jj)+uij (ii-1, jj)) 
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+ (uij(ii,jj+l)-2*uij(ii,jj)+uij(ii,-jj-l))) 
+alpha*dx~2*(uij (ii, j j)-uij (ii, j j)~3) ; 
°/„+eps*dx~2*rr*uij (ii, j j) ; 

°/„+eps*dx"2*sub({xi=(ii-n-l)/n,yi=(j j-n-l)/n},f f ij) ; 

o o o 

Last is the amplitude condition that the central value is the unchanged 
amplitude (could make the amplitude equal to the element mean if neces- 
sary), and set ok if all the right-hand sides are zero. 

» ncsm2d «+ 

eqns(q:=q+l,l) : =ud(n+l ,n+l) ; 
ok:=if eqns=zeros then 1 else 0; 

o o o 



Form matrix of the linear equations In the first iteration only, now 
form matrix of the linear equations from the perturbation variables ud and 
gd that have been carried around in this first iterate. Extract these unknowns 
and assign their coefficients to the lu matrix. Then perform the LU decom- 
position for later solution of the homological equation. Lastly remove the 
symbolic unknowns. 

» ncsm2d «+ 

if iter=l then begin 

for qq:=l:nn do begin q:=0; 

for ii:=l:ns do for jj:=l:ns do 

lu(qq,q:=q+l) :=-coef fn(eqns(qq, 1) ,ud(ii,jj),l); 

lu(qq,q:=q+l) : =-coef f n(eqns (qq, 1) ,gd, 1) ; 

end; 

write "starting LU decomposition to store"; 
lu_decomp() ; 

write "finished LU decomposition"; 
showtime ; 

Comment get rid of perturbation variables hereafter; 
let { ud(~ii,~jj)=>0, gd=>0 }; 
end; °/> the first iterate linear equation stuff 

o o o 

Use the factorisation stored in lu to solve for the updates, and assign to 
the relevant g and Uy. 

» ncsm2d «+ 
write "Back substitution solves for updates"; 
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lu_backsub() ; 

gij :=gij+eqns(nn, l)/dx"2; 

q:=0; 

for ii:=l:ns do for jj:=l:ns do 

uij (ii.jj) :=uij (ii , j j )+eqns(q: =q+l , 1) ; 

ooo 

Terminate the loop when all residuals are zero. 

» ncsm2d «+ 

showtime ; 

if ok then write iter : =1000000+iter ; 

end; 

ooo 



5 Postprocessing numerical 

5.1 Write a sci/matlab function script 

Optionally generate an efficient sci/matlab function in file gl2dnumred. This 
function file is to be processed with the unix script red2mat and then it can 
be used in matlab to integrate solutions using ode solvers. 

» ncsm2d «+ 

write " 

Generating mat/scilab function file in gl2dnumred 
ii . 

« functionfile » 

write "finished generating mat/scilab file"; 

write "Use sed -f red2mat gl2dnumred > gl2dnpde.m"; 

showtime ; 

ooo 



5.2 Derive the equivalent PDE 

Optionally derive the equivalent PDE for this model. 

» ncsm2d «+ 
if 1 then begin write " 

Optionally find the equivalent partial differential equation 
of the holistic discretisation. 
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« equivpde » 
showtime ; 

end; 

ooo 



5.3 Derive the centred difference form 

Optionally derive the finite difference form of the evolution. 

» ncsm2d «+ 
if 1 then begin write " 

Generate evolution in terms of centred finite difference 
operators md=mu*delta and dd=delta~2 in each of the x and 
directions . 

« finitediff » 
showtime ; 

end; 

ooo 

Additionally write the coefficients of three parts in the C(tx 3 +y 3 ) model. 

» finitediff «+ 

write g2nd:=f actorize(coeffn(h~2*gop,gam,2)) ; 
write g3rd:=f actorize(coeffn(h~2*gop,gam,3)) ; 
write gnon:=f actorize(coeffn(coeffn(h~2*gop,alf , 1) ,gam, 1)) 

ooo 

Gather results of gnon into a list and fit rational Pade functions in 1 /n 2 . 
Need an odd number of results. Note that rn = 1/n and the multiplication 
by four of the coefficient. 

» finitediff «+ 
write "Some meta-analysis of different subgrid resolutions 
ns:={2,3,4,5,6,7,8}; 
cs:={l/72 
,1/60 

,179/10136 

,775/42762 

,679909/36998632 

, 237808723/12834019900 

, 133046058951/7141880630840 

}; 
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operator a; operator b; 
procedure paden(ns , cs) ; 
begin 

o:=(length(ns)-l)/2; 

vars:=for p:=0:o join {a(p) ,b(p)}; 

eqns : =(b(0)=l) . (f or 1 : =1 : length(ns) collect 

(for p:=0:o sum b(p)/part (ns , 1) ~ (2*p) )*part (cs , 1)*4= 
(for p:=0:o sum a(p)/part (ns , 1) " (2*p) ) ); 
soln : =solve (eqns , vars) ; 

return sub(soln, (f or p:=0:o sum a(p)*rn~ (2*p) ) 

/(for p:=0:o sum b(p)*rn~(2*p))) ; 

end; 

padefull:=paden(ns,cs) ; 

on rounded; write largenlimit :=sub(rn=0,padefull) ; off rounded; 
write pade4:=paden({2,3,4}-, 

{1/72 

,1/60 

,179/10136}) ; 

o o o 

5.4 Fin numerical 

Give the grand final 'end' statement, and output all the code. 

» ncsm2d «+ 

end; 

o o o 



6 Factorisation and solution of linear equa- 
tions 

6.1 LU decomposition 

This LU decomposition was adapted from Numerical recipes in Fortran 77 [9]. 
Requires matrices lu(n,n), indx(n, 1) and eqns(n, 1) to be predefined. Ma- 
trices are global in scope in Reduce so we do not bother passing parameters. 

» lu_decomp « 

% see cadsmpd2d.pdf for documentation, adapted from Press et al . 
procedure lu_decomp; 
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begin scalar n,np,i, j ,k,imax,dd, aamax , dum , sum, tiny; 
tiny:=1.0e-20; 
n:=f irst (length(lu) ) ; 
dd:=l; 

for i:=l:n do begin 
aamax :=abs (lu(i , 1) ) ; 

for j:=2:n do aamax : =max (aamax, abs(lu ( i, j))) ; 
vv(i , 1) :=l/aamax; 
end; 

for j:=l:n do begin 
for i :=1 : (j-1) do 

lu(i, j) :=lu(i, j)-(for k:=l:(i-l) sum lu(i,k)*lu(k, j)) ; 
aamax : =0 ; 

for i:=j:n do begin 

lu(i, j) :=lu(i, j)-(for k:=l:(j-l) sum lu(i,k)*lu(k, j)) ; 
dum: =vv(i , l)*abs(lu(i , j ) ) ; 
if dum>=aamax then begin 
imax:=i ; 
aamax : =dum ; 
end; 
end; 

if not(j=imax) then begin 
for k:=l:n do begin 

dum : =lu ( imax , k) ; 

lu(imax,k) :=lu(j ,k) ; 

lu( j ,k) : =dum; 
end; 

dd:=-dd; 

vv(imax, 1) :=vv(j , 1) ; 
end; 

indx(j , 1) :=imax; 

if lu(j,j)=0 then lu(j , j) :=tiny; 
if not(j=n) then begin 
dum:=l/lu(j , j) ; 

for i:=(j+l):n do lu(i, j) :=lu(i, j)*dum; 
end; 
end; 
end; 
end; 

o o o 
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6.2 LU back substitution 



This LU back substitution was adapted from Numerical recipes in Fortran 77 [9]. 
Requires matrices lu(n,n), indx(n, 1) and eqns(n, 1) to be predefined: the 
first two must be obtained from procedure lu_decomp; and the last may be 
algebraic. Matrices are global in scope in Reduce so we do not bother passing 
parameters. 

» lu_backsub « 

°/„ see cadsmpd2d.pdf for documentation, adapted from Press et 

procedure lu_backsub; 

begin scalar n,i, j ,ii,ll,sum; 

n:=f irst (length (lu)) ; 

ii:=0; 

for i:=l:n do begin 

ll:=indx(i,l) ; 

sm: =eqns (11 , 1) ; 

eqns (11 , 1) : =eqns (i , 1) ; 

if not(ii=0) then 

sm:=sm-(for j:=ii:(i-l) sum lu(i , j ) *eqns ( j , 1) ) 

else if not(sm=0) then ii:=i; 

eqns (i , 1) :=sm; 
end; 

for i:=n step -1 until 1 do eqns (i , 1) : =(eqns (i , 1) 
-(for j:=(i+l):n sum lu(i , j )*eqns (j , 1) ) )/lu(i , i) ; 

end; 
end; 

o o o 
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